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Regularized Non-Gaussian Image Denoising 
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Abstract 

In image denoising problems, one widely-adopted approach is to minimize a regularized data-fit ob¬ 
jective function, where the data-fit term is derived from a physical image acquisition model. Typically 
the regularizer is selected with two goals in mind: (a) to accurately reflect image structure, such as 
smoothness or sparsity, and (b) to ensure that the resulting optimization problem is convex and can 
be solved efficiently. The space of such regularizers in Gaussian noise settings is well studied; however, 
non-Gaussian noise models lead to data-fit expressions for which entirely different families of regularizers 
may be effective. These regularizers have received less attention in the literature because they yield 
non-convex optimization problems in Gaussian noise settings. This paper describes such regularizers and 
a simple reparameterization approach that allows image reconstruction to be accomplished using efficient 
convex optimization tools. The proposed approach essentially modifies the objective function to facilitate 
taking advantage of tools such as proximal denoising routines. We present examples of imaging denoising 
under exponential family (Bernoulli and Poisson) and multiplicative noise. 

Keywords. Image denoising, convex optimization, total variation, sparsity, proximal operators, image 
regularization 


1 Introduction 

In image reconstruction, measured image data is always corrupted by noise. Active research addresses 
the problem of removing various types of noise from the observations based on models of both the data 
acquisition process and of underlying image structure. In many settings this noise is modeled using a Gaussian 
distribution because it is a good model of errors associated with common detectors like CCD arrays. However, 
many more specialized applications yield noise better modeled with alternative distributions, such as Poisson 
HU, Gamma m, Rayleigh m, Bernoulli [20], and Geometric m- In these and other settings, a common 
denoising approach is to find an image which minimizes a regularized negative log-likelihood of the noisy 
observations, with regularizers adopted from the Gaussian image denoising literature. This paper shows 
that alternative families of regularizers, adapted to the noise model, can yield empirical performance gains. 
Furthermore, regularizers which are non-convex in Gaussian settings can be convex and quite effective in 
these non-Gaussian settings. 

More specifically, our interest is in recovering a true underlying image x* G from the noise-corrupted 
observations y G where y is drawn from a distribution parameterized by x* which models the noise, 
with the density denoted p{y\x). To obtain an estimate x G of x*, a common approach is to minimize a 
smooth, convex data-fit term (/>(•) regularized by a convex (and potentially non-smooth) regularization term 
/-(O To formalize this notion, we consider 

X = arg min {ffi(x, y) := 0(x, y) + Tp(x)}, 

a: 

where r > 0 is a tuning parameter. In Gaussian noise settings, we often use (t){x^y) = \\x — yWl^ but more 
generally, a common choice for (t){x^y) is a term proportional to — logp(^|x), the negative log-likelihood of y 
given X according to a distribution selected to model noise statistics. In general, (j) reflects the characteristics 
of the noise corrupting our observations. 

Wniversity of Wisconsin-Madison, Department of Electrical and Computer Engineering 
1415 Engineering Drive, Madison, WI 53706. 

^Kernel and patch based methods, including bilateral filters, non-local means, and BM3D implicitly employ a data-dependent 
regularization function. Such approaches are beyond the scope of this paper. 
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The regularize! p(') has been extensively studied and there are a few examples widely used in image 
recovery techniques. The choice of this regularize! depends on the assumptions made about the underlying 
image structure. Popular choices include the £i norm for coefficient sparsity in a wavelet basis or other 
dictionary, the total variation (TV) semi-norm for image gradient sparsity, and Huber-like functions which 
are akin to the ii norm but smooth. In general, image processors choose a regularize! according to two 
desiderata: 

(Dl) the objective function ^{x^y) may be minimized efficiently; and 
(D2) the regularize! p{x) accurately reflects image structure. 

For instance, when the data-fit term is simply the sum of squared errors — ^||25 fhen one would choose a 
convex regularize! p. In this paper, we describe several classes of image denoising problems where 
the noise is non-Gaussian and where non-convex regularization functions (which would be 
undesirable in Gaussian noise settings) not only admit efficient convex optimization algorithms, 
but also yield substantial empirical performance gains. 

As an illustrating example, consider two different denoising approaches, both based on total variation 
regularization mi: 


ii=argmin -logp(2/|a:) + n||a:||TV 

X 


and 


X2 =arginin - logp(y|a;) + T 2 II loga;||Tv, 

X 


( 1 ) 

( 2 ) 


where 


Pl-lp2-l _ 

\\x\\tv '-E E y + (Xij - Xij+i^ 

i=l j=l 

Pi-1 P 2 -I 

^ ~ ^i-\-i,P2\ V \^pid ~ ^pi,j+i\ (3) 

i=i j=i 

is known as the discrete isotropic TV [2] , where in a slight abuse of notation we consider x to be an pi x p 2 

image matrix. 

Comparing the two optimizations, the main difference is that the TV norm is applied to either the image 
[8] or the log of the image (assuming all the pixel values are positive) [TT]. If we think about the TV norm 
as promoting sparse image gradients, then it is clear that applying the TV norm to the log-intensity image 
is a sensible reflection of image structure, and may allow for larger non-zero gradient magnitudes than the 
TV norm applied directly to the intensity. If the noise were Gaussian, so that — logp{y\x) oc ||^ — 
the optimization problem in Q would be non-convex, making this regularize! undesirable. However, in 
the Poisson noise setting, where — \ogp{y\x) oc + Vi ^ straightforward reparameterization (as 

described in m) leads to a convex optimization problem. 

In other words, the shape of the log-likelihood funetion helps determine the set of regularization funetions 
that satisfy (Dl) and (D2). This means that the best ehoiee of regularizer for a given problem may depend on 
the noise model and henee may be quite different from what one would ehoose in a Gaussian noise setting. 


1.1 Organization of the paper 

Section contains an explicit formulation of the the problem, and Section details our approach. Section 
presents how reparameterization can affect computational-statistical tradeoffs. In Section we introduce a 
generalization for image reparameterization for the exponential family of distributions; in Sections 5.1 and 
|5.2| we consider two specific examples of Bernoulli and Poisson noise, respectively. In Sectionwe consider 
an example of reparameterization for multiplicative (speckle) noise. In Section]^ we present experimental 
results in selecting the best reparameterization (in terms of RMSE) over a bank of images. We conclude the 
paper in Section 
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2 Problem formulation 

In choosing a good regularizer, one might search over a space of candidate regularizers such that the ob¬ 
jective function ^{x,y) = (j){x^y) -h p{x) is convex and choose a regularizer p from this space which yields 
accurate estimates across a collection of images. However, this space can be very difficult to specify and 
overly restrictive in practice, especially since (t){x^y) is not necessarily convex in x. We consider a tractable 
alternative in this paper. Specifically, we start with a fixed^ known family of convex image regularizers p, 
denoted IZ. For each p e we search over a space of candidate reparameterization functions f : W W 
such that the objective function 

^(/(^')> y) = + Tp{0) 

is convex in 0 for all r > 0. For any such /, given a new parameter 0 ^ x = f{0) G is a reparameterized 

version of the image intensity. In other words, instead of searching over regularizers that satisfy (Dl) and 
(D2), we search over corresponding reparameterizations of the log likelihood. This approach guarantees that 
we will be able to leverage off-the-shelf convex optimization solvers that rely on proximal operators for a 
regularization function p 

More specifically, let 7^ be a collection of convex regularization functions p :W^ that reflect desirable 
image structure {e.g., sparsity, smoothness, etc.). Let 0 C be a feasible parameter space. Define the set 
of reparameterization functions 

= {/ : 0 : 0(/(-),p) convex and smooth}. 


Given some / G p G 7^, and r > 0, we may estimate our image via: 


0f,p,T = argmin {i>(/(6»), y) = (p{f{e),y) + Tp{e)} 
9 


( 4 ) 


where it is implied that and Xf^p^^ are functions of the noise p, but we omit writing this explicitly 

throughout this paper for brevity. 

Our goal is to choose a regularizer p G 7^ and reparameterization f G J- according to 

{r^pl= argmin R{f,p) 

{f,p)eTxn 

where R measures the risk 


R{f.p) 




min 

r>0 




^I|2/||X||2 


( 5 ) 


where the expectation is with respect to both the noise in the data y (of which Xf^p^r is an implicit function) 
and a distribution over the space of images. 

Notice that because the variable we are optimizing over is we are directly regularizing the new parameter 
0 and not the image. The interpretation is that by changing the reparameterization function /, we change the 
relationship between the image x and the parameter vector and thus implicitly change how x is regularized. 
To see this directly, note that if / is invertible, is equivalent to 

Xf,p,r = argmin {^{x,y) = (j>{x,y) + Tp{f~'^{x))}. (6) 

X 

We note that it is not required for / to be invertible in general, but it is useful in establishing this inter¬ 
pretation. From this perspective, it is clear that the regularizer applied to the image is p(/“^(x)), and so 
choosing a reparameterization function / implicitly selects an associated regularizer. Furthermore, the im¬ 
plicit regularization function p(/“^(')) is not necessarily convex, but convex optimization methods can still 
be used to denoise the image. 

The key point is that the form of the log-likelihood will govern the eontents of R and the speeifie risk¬ 
minimizing reparameterization and assoeiated regularizer. Conventional wisdom tells us that a “good” reg¬ 
ularizer should simply be convex and reflect desirable image structure, and once a good regularizer is found 
it can be used with any noise model. In contrast, we will see examples where regularizers considered good 
according to this conventional wisdom are suboptimal. Furthermore, regularizers identified in our framework 
will often be non-convex functions of the image (and hence undesirable in Gaussian noise settings) but still 
leading to convex optimization formulations for non-Gaussian image estimation. 
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3 Proposed approach 

Our goal is to minimize the risk defined in (§. Since in general we do not have access to the distribution 
over the space of images, we minimize an empirical estimate of risk using a collection of n representative 
images, denoted ..., For the image, we generate m noisy realizations ..., whose 

reconstructions are denoted • • •, via Depending on the noise distribution assumptions, we 

select an appropriate data-fit term (j) based on the negative log-likelihood. Note that often the negative log- 
likelihood itself will be convex; however, alternative formulations of 0 may be desired when it is non-convex; 
we present an example of this in Section]^ For a particular (/, p) pair, we then compute the empirical risk 


R{f,p) 


n 


m 


nm 


E““Ei 

i=i j=i 




r*(d| 


(7) 


which is essentially a metric that averages the normalized RMSE of n image examples with m noise realizations 
each, noting that r is chosen to minimize the average (across noise realizations) RMSE for each image. Einally, 
we are then interested in choosing a reparameterization / that minimizes this empirical risk, namely 

/ = argmin R{f). 
feT 


In order to compute the estimates we adopt a proximal gradient method m, where a key step in 

iteration t is to solve the proximal denoising problem 


= argmin i||6> - St\\l + ^p{0) 


at 


(8) 


for appropriately specified gradient descent step St and (potentially adaptive) step size at. Because our 
strategy is to search over reparameterization functions, off-the-shelf proximal gradient methods can be used 
to compute quickly and accurately P[6l 1^ [22] . 


4 Computational-statistical tradeoffs 

The proposed methodology can be used not only to improve empirical RMSE performance, but also to exploit 
computational-statistical tradeoffs main]. Take the example of Poisson denoising in Section comparing 
0 and 0. In light of Section]^ we can rewrite 0 equivalently as 

62 = argmin - logp{y\ exp(6»)) + T2||6 »||tv 
9 

X2 = exp (^2), 

where exp(-) is applied element-wise, which is how we would perform the “logarithmic regularization” of 0 
in practice: first reparameterize the image with an exponential transform to obtain the estimate 82 ^ then 
obtain the image estimate through the reparameterization. With this optimization problem, the data-fit 
term — logp(^| exp6>) oc exp6>^ — ytOi grows without bound (i.e. does not have a Lipschitz gradient) as 

0 ^ 00 but is otherwise well-behaved (i.e. does have a Lipschitz gradient). On the other hand, — log p{y\x) oc 
— yilogXi grows without bound as x ^ 0 but is otherwise well-behaved. So we see that in both 
cases, we have data-fit terms whose gradients are not globally Lipschitz, but are locally Lipschitz on some 
domain. 

A Lipschitz gradient is important in characterizing the convergence rate of many convex optimization 
methods [3]. In general, we can say that an algorithm applied to a function with a smaller Lipschitz constant 
on its gradient will converge much faster than one with a larger Lipschitz constant; as a general intuition, 
having a larger Lipschitz constant (i.e. greater curvature) requires smaller steps in a gradient descent setting, 
leading to slower convergence. This leads to a practical problem: many of the data-fit terms resulting from 
commonly-used noise distributions, and not just Poisson, often have regions with unbounded curvature. 
The class of reparameterizations R described above can be constrained to yield data-fit terms with globally 
Lipschitz gradients, which in turn increases convergence speed; in this case f{0) is a function that is selected 
to make regions of (t){f{0)^y) with unbounded curvature in f{0) have bounded curvature in 0 . 
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In settings where convergence speed is critical, we can use this idea to modify the objective function 
used to select a reparameterization and reguarization pair. Specifically, let L{(j){f^y)) denote the Lipschitz 
constant of the gradient of Then instead of minimizing R{f,p) as in 0 , we might compute 

{f*,P*)= argmin R{f,p) + X{L{(p{f,y)))~^, 

(f,p)eJ^xn 

where A > 0 controls the computational-statistical tradeoff. This idea will be demonstrated in the results 
presented in Section 


5 Reparameterization for the exponential family 


Noise from the exponential family is commonly seen in imaging problems, so it is useful to present a general 
0(/(^),^) function for this class of problems. Distributions in the exponential family can be written in the 
form 

p 

p{y\v) = n exp [riiT{yi) - A{rii)], 

i=l 


where h(') is the base measure (a constant to scale the distribution so that it integrates to 1), y e W are 
the natural parameters that depend on distribution parameters, T(') is the sufficient statistic, and A{') is 
the log-partition function [T]. The negative log-likelihood of this family is then proportional to 

p 

Hv, y) = 'E “ ViT{yi)- (9) 

i=l 

Equivalently, we can rewrite the data-fit term 0 in terms of the image intensities x, whose relationship to 
the natural parameter is y = g{x), where ^ i-^ is an invertible function that maps from the image 

X to the natural parameter 77 , defined for all distributions in the exponential family; we note that g cannot 
always be applied element-wise, but in the specific examples to follow, we use distributions whose g functions 
can be applied element-wise. We may then write the negative log-likelihood as 


(Pix,y) = '^A{[g{x)]i) - [g{x)]iT{yi). 


( 10 ) 


To develop our family of reparameterizations we consider functions which bridge between the data fit 
terms in ^ and ( 10 ) via 


Hf{o),y) = E^wm]i) - wmmyi). 


( 11 ) 




Our task now is to form a class of reparameterization functions E where the following criteria are met: 

1. for each f ^ E, the corresponding 0(/('),^) is continuous, convex, and has Lipschitz gradients. 


2 . regularizing in the new parameter space yields strong empirical performance; and 

3. we may efficiently explore the space E. 

To accomplish this, given a, 6 G M, define 

fa,h(P) = g~^[a^ + h). 

We define a class E where [f{0)]i is either [fa,b{^)]i or depending on which has a smaller partial derivative 
with respect to Oi. Formally, 


E = 


f continuous over 0 : 3a, 6 G M s.t. 


[f{0)h = 


[faA0)h 

0^, 


d[fa,bmi 

dOi 


< 1 


otherwise 



( 12 ) 
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Note that for / G 


[gifiom 


dOi 6 , 

[9{S)]u 


9[fa.bmi 

dOi 


< 1 


otherwise 


(13) 


We note that this f{0) reparameterization does not always yield a convex data-fit term for all distributions 
in the exponential family, though in what follows, we explore examples that are convex. We also note that 
the family defined in is only one possible family and others are certainly possible, though we found 
this family to be useful for the two examples in the following subsections. In both of the examples to 
follow, it is straightforward to show that the reparameterized data-fit terms we propose are convex and 
have Lipschitz gradients through basic algebra (convexity by showing that the diagonal Hessian matrix has 
diagonal elements ^ nonnegative, and Lipschitz gradients by showing are finite); 

continuity is guaranteed by construction. 

Table [2 lists the relevant functions for the Bernoulli and Poisson distributions, which will be discussed in 
the following. 



A{-) 

T{y) 

g{-) 


Bernoulli 

log(l + exp(-)) 

y 

log(i3:) 

exp(-) 

l+exp(-) 

Poisson 

exp(-) 

y 

log(-) 

exp(-) 


Table 1: Table of relevant functions for Bernoulli and Poisson distributions (log-partition H, sufficient statistic 
T, natural parameter mapping g, and inverse mapping g~^). In these specific examples chosen for this paper, 
all the operations in this table can be applied element-wise. 


5.1 Bernoulli noise (noisy one-bit sensors) 

One example of a denoising problem that fits into this framework is in quantization noise removal. Consider 
an example of an unreliable sensor system where measurements are known to be imprecise, and so the 
observations are taken to be binary m- In this one-bit quantization noise case, the sensor can be modeled 
with a dithering and quantization so that 


Vi = C 


Xi >Ui 


1 , if Xi > Ui 
0 , otherwise 


where Xi G [0,1] and corresponds to relative pixel intensities, and Ui is a threshold drawn from Uniform(0,l). 
This can be modeled with a Bernoulli distribution, namely y ^ Bernoulli(x*), which falls into the exponential 
family. 

Adopting the strategy above in ( 11 ) and applying Tablethe Bernoulli distribution results in the data-fit 
term 


= E + exTp[g{f{0))]i) - yi[g{f{0))]i. 


i=l 


Before defining f{0) for the Bernoulli case, we note that we use different reparameterizations depending on 
if Oi < 0.5 or Oi > 0.5; this is a small extension of (12). For this particular case we specify a family of 
reparameterizations as: 


E = < f continuous over 0 : 3a, 6 i, 62 G M s.t. 


[fmi = 


' [faM (^)]i < 1 and Si < 0.5 

< 1 and > 0.5 , i = 1, • • • ,P 

Oi otherwise 
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Note that for any function in this class, g{f{0)) has the form: 


wmi = < 


adi + h, 


< 1 and Oi < 0.5 


adi + b- 2 , < 1 and 9, > 0.5 ■ 


'^2> 00. 
otherwise 

In the above, there are three degrees of freedom in a, 6 i, and 62 - However, there is effectively only one 
because for any triple (a, 5i, 62 ) satisfying the continuity constraints, there exists some k G (0, 0.5) such that 


[gifiom = 


fe(l-fe)^* lb’ ^ ^ 

log(yz^), k<ei<l-k , 

+ ei>i-k 

where image intensities are remapped via 


if mi ■■= 


^ —h 

Vexp 


(/c(l-/c)^* + l-/c) 


1 -1 


Oi, 


-1 


l + ^exp(^^^, + i) 
so that the family of reparameterizations can be formally rewritten as 


6 i < k 

k < Si <l-k , 
0 i>l-k 


(14) 


< 


/:e 


: 3k e (0,0.5) s.t. 



1 + exp 1 

^/c(l-/c)^* + 1 -k 

:)] 

-1 

5 

0 i<k 


[Oi, 




k < 9i <1 — k , i = 1,- ■ 



t + T^exp| 

^/c(l-/c)^* k) 


1 

5 

9i> 1-k 

> 


Summarizing, in the Bernoulli case we propose the data-fit expression 




(15) 




where 

4>i{fW,y) ■= ' 

for some k G (0, 0.5) 


'log (1 + exp(^^)) - Vi 

log(T^) -2/ilog(Y^), 

log (1 + exp(|^)) - 2/i 


fe(i-fc) +log(T^) 

'l^^+log(^) 


Oi < k 

k < Oi < 1 - k 
0 i>l-k 


It is straightforward to show that the reparameterized data-fit term (j){f{0)^y) we propose is convex and 
has Lipschitz gradients through algebra (convexity by showing that the diagonal Hessian matrix has diagonal 
elements ^ nonnegative, and Lipschitz gradients by showing are finite). 

Our proposed f(0) is an invertible function. Invertibility is not required in our setting, but it is useful in 
understanding the implicit regularization function p{f~^{x)) seen in (§. The inverse function is 

-log(T^(^-l))] fc(l-fc), Xi<k 
{f~^ix)]i = {xi, k < Xi < 1 - k . 

i - log (^(^ - 1)) k{l-k), Xi>l-k 


( 16 ) 
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0 0.25 0.5 0.75 1 

Pixel value 1 



0 0.25 0.5 0.75 1 

Pixel value 1 



0 0.25 0.5 0.75 1 

Pixel value 1 


(a) Intensity reparam. 


(b) Natural reparam. 


(c) Proposed reparam. (k = 0.4) 


Figure 1: Bernoulli implicit regularizers: In these illustrations, each axis represents varying pixel values 
and the brightness shows the absolute difference |[/~^(^)]i — [/~^(^)] 2 | for: 


(a) 

(b) 

(c) 


/ ^{x) = X (using the f{0) = 9 pixel intensity reparameterization) 


/ ^{x) = log(j^) (using the /(0) =g ^{9) = 


_ exp 9 

1+exp 9 


natural reparameterization) 


/ ^{x) as defined in (IM for k = 0.4 (using the proposed reparameterization (14)). 


Because the TV norm essentially takes image gradient differences, these illustrations give an intuition behind 
the implicit regularizers introduced by each reparameterization when using p(') = || • ||tv- We see that when 
pixel values are within [0, k] and (1 — k, 1] our proposed regularizer in (c) behaves like (b), but for pixel values 
in (/c, 1 — k] our proposed regularizer (c) behaves like (a); this is consistent with our piecewise construction 
of our proposed reparameterization. 


For an intuition behind the regularization. Figure illustrates how the implicit regularizer ||/ ^(x)||tv 
behaves. 


5.2 Poisson noise 

Poisson noise is a problem that is of wide interest in the image processing community, where it is also 
sometimes referred to as shot noise. This noise is prevalent in low-light imaging where the dominant noise 
source comes from photon counting; this arises in various fields of study where images are acquired in dark 
backgrounds (i.e. night vision) or are acquired with a short exposure time (i.e. biomedical imaging). In 
these photon-limited realms, the noise is modeled as being drawn from a Poisson distribution whose 

rate parameter represents the true underlying signal intensity, namely y ^ Poisson(x*). 

In the context of our generalized model from (11), applying Table leaves us with 




i=l 


One candidate g{f{0)) function comes from the generalization (13). For E as defined in (12), 


[»(/(«))!.+ 

[log(6>i), 9i> ^ 

where the image intensities are mapped back for some k G (0, oc) via 

f expikOi + log(b) - 1), 9i< I 


imh ■■= 


I 9i, 




which is similar to the function proposed in [12]. In applying (12) to the Poisson case, there are two degrees 


of freedom in a and h. However, there is only one degree of freedom in the above expressions as a result of 
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the pair (a, b) satisfying continuity constraints. The family for reparameterizations can be formally rewritten 
as 


T = 


f:e^W:3k>0 s.t. [f{9)]i 


exp(A:6'i + log(i) - 1), Oi < I 

^*5 k 


i = I,--- ,p 


5.2.1 Alternative Poisson reparameterization 

Another candidate reparameterization to explore is based on a function already proposed in the literature in 
[16], where 


9 {fi^)) ■= log(fclog(l + exp(6»))), 

where the image intensities are remapped via 

f{0) := k log(l + exp(0)), (17) 

in what we call the “log-exp” parameterization. Though this function was motivated by a Poisson matrix 
completion problem, it was designed in order to address the unbounded curvature of the data-fit when 
estimating the rate parameter directly. We introduce the k parameter in this work to describe a family of 
reparameterizations, which can be formally stated as 

E={f-.e^W : 3k >0 s.t. f{d) ■- fclog(l + exp(6»))}. 


Summarizing, in the Poisson case we propose the data-fit expression 

p 

^ + exp(0i)) - Vi log(fc log(l + exp(6»j))) (18) 

i=l 

for some k G (0, oc). 


It is straightforward to show that the reparameterized data-fit term (j){f{0)^y) we propose is convex and 
has Lipschitz gradients through algebra (convexity by showing that the diagonal Hessian matrix has diagonal 
elements ^ f^at are nonnegative, and Lipschitz gradients by showing are finite). 

Our proposed f(0) is an invertible function. Invertibility is not required in our setting, but it is useful in 
understanding the implicit regularization function p{f~^{x)) seen in (§. The inverse function is 


/ ^(x) = log(exp(-) - 1). 


(19) 


For an intuition behind the regularization, Figure illustrates how the implicit regularizer ||/ ^(x)||tv 
behaves. 


6 Reparameterization for speckle noise 

Speckle noise is found in many imaging applications, particularly in Synthetic Aperture Radar (SAR) imag- 
ing [13] and ultrasound [9]. Speckle differs from other noise types in that it is multiplicative, where the 
observations obtained from an imaging system are the products of the true underlying intensity and draws 
from a Gamma distribution. The observation model can be written element-wise as follows: 

Vi ^ 2^2 5 

where n is multiplicative noise that is known to be from a Gamma distribution with density 

Virii) = |;^<“^exp(-/3ni), 
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Pixel value 1 



0 12 3 4 

Pixel value 1 



Pixel value 1 


(a) Intensity reparam. 


(b) Natural reparam. 


(c) Proposed log-exp (k = 2) 


Figure 2: Poisson implicit regularizers: In these illustrations, each axis represents varying pixel values 
and the brightness shows the absolute difference |[/~^(^)]i — [/~^(^)] 2 | for: 

(a) f~^{x) = X (using the f{0) = 0 pixel intensity reparameterization) 

(b) f~^{x) = log(x) (using the f{0) = g~^{0) = exp(^) natural reparameterization) 

(c) f~^{x) as defined in for /c = 2 (using the proposed reparameterization 0)- 


Because the TV norm essentially takes image gradient differences, these illustrations give an intuition behind 
the implicit regularizers introduced by each reparameterization when using p(-) = || • ||tv- We see that when 
pixel values are small the proposed regularizer in (c) behaves like (b), but for pixel values that are large our 
proposed regularizer (c) behaves like (a). 


where a and P are shape and rate parameters, respectively. In SAR imaging, for example, it is commonly 
assumed that a = ^ = S, where S is the number of “looks” (or measurements) [19]. Using this, we find that 
the negative log-likelihood p{y\x) is then proportional to 

4>{x,y) = y^S'logXj + S—. 

i=i 

Reparameterization results in the data-fit term 

Hfm,y) = p^S\og[f{9)]i + 

We first consider the reparameterization 

[hid)]i := l/Oi 

that yields the data-fit ~ fog^i+^i^i? which is convex; this function, however, has unbounded curvature 

as ^ 0. To remedy this, we can use an exponential transform via 


[f2{d)]i ■■= 1/exp(6»i) 

to obtain —9i + y^exp^j; this is also convex, but has unbounded curvature as 0 —> oo. To address the 

unboundness of both possible solutions, we combine them through 


umi 


k exp(0i) ’ 

. fc(ibi)’ 


ei<o 

0i>0' 


( 20 ) 


where k > 0 affects the curvature, and allows us to define a family of reparameterizations over which to 
search. Formally, the family this describes is 




3k > 0 s.t. [f{0)]i := 


-{ 


k exp(0i) ’ 
1 

k{l-\-6i) ’ 


0 i<0 

6»i > 0 ’ 
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0 12 3 4 

Pixel value 1 



0 12 3 4 

Pixel value 1 



Pixel value 1 


(a) fi{0) reparam. 


(b) 12 ( 0 ) reparam. 


(c) Proposed reparam. (k = 2) 


Figure 3: Speckle noise implicit regularizers: In these illustrations, each axis represents varying pixel 
values and the brightness shows the absolute difference | [/' - [/ ^(^)]2| for: 


(a) [f~^{x)]i = ^ (using the [fi{0)]i = ^ reparameterization) 

(b) = — log(x) (using the [f 2 { 0 )]i = reparameterization) 

(c) as defined in ( [2^ for /c = 2 (using the proposed reparameterization ([2Q|)). 


Because the TV norm essentially takes image gradient differences, these illustrations give an intuition behind 
the implicit regularizers introduced by each reparameterization when using p(-) = || • ||tv- We see that for 
small pixel values within [0, the proposed regularizer in (c) behaves like (b), but for large pixel values within 
(^,+oc) our proposed regularizer (c) behaves like (a); this is consistent with the piecewise construction of 
our proposed reparameterization. 


Summarizing, in the speckle noise case we propose 






( 21 ) 


where 




yikexp{9i) - 9i - log k, < 0 

yik{l + 9i) - log(l + 9i) - log k, > 0 


for some k G (0, oc). 


It is straightforward to show that the reparameterized data-fit term (j){f{0)^y) we propose is convex and 
has Lipschitz gradients through algebra (convexity by showing that the diagonal Hessian matrix has diagonal 

elements ^ f^at are nonnegative, and Lipschitz gradients by showing are finite). 

Our proposed f{0) is an invertible function. Invertibility is not required in our setting, but it is useful in 
understanding the implicit regularization function p{f~^{x)) seen in §. The inverse function is 




L kxi 



( 22 ) 


For an intuition behind the regularization, Figure illustrates how the implicit regularizer ||/ ^(a:;)||TV 
behaves. 
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7 Numerical experiments 

As described in Section we seek the reparameterization / whose resulting estimate x minimizes the average 
normalized RMSE (||x — x* || 2 /||x* II 2 ) over a family for a bank of images for a given noise level. In other 
words, we are interested in finding the reparameterization / that minimizes the empirical risk R(/, p) as in 
0 , or 

f = avg min R{f,p). (23) 

feT 


To be clear about the experimental results to follow (and to tersely sum up our discussion in previous 
sections), we obtain the estimatesfor the image and noise realization by using: data-fit (15) for 
Bernoulli noise, log-exp data-fit ( p^ for Poisson noise, and data-fit (21) for speckle noise. To showcase the 
value of our reparameterizations, we compare our experimental results to a baseline denoising that minimizes 
the negative log-likelihood with a regularizer applied directly on the image intensities, namely. 


Xr = argmin {^{x,y) = (l){x,y) +rp(x)}, (24) 

X 

which is parameterized by the tuning weight r. 

In the examples to follow, we use the TV regularizer as defined in 0 for p. We note that all r weighting 
values for the reparameterization and baseline estimates in the following experiments were chosen clairvoy- 
antly, as to minimize each realization’s RMSE. 

Eor our proximal gradient method, we adopt SpaRSA m to compute the minimizer of -h 

t|| 6 >||tv 5 foi* fhe image and its noise realization. The main idea of the method is to iteratively perform 
0 , where St = 0^ — ^V(/)(/( 6 >^), is a step in the direction of steepest descent and at is a step-size that 

changes at each iteration. Eor a TV regularizer, ^ can be solved with the East Gradient Projection (EGP) 
algorithm [ 2 ]. We defer more detailed descriptions of both SpaRSA and EGP to their original papers, as 
they are not the main focus of this work but are used for experimental validation. 


7.1 Shepp-Logan phantom 

As a first example to demonstrate the value of reparameterizing, we show example reconstructions with the 
Shepp-Logan phantom (512 x 512) for each of the noise cases considered. We note that for Bernoulli noise 
ground truth was scaled to have intensities between 0 and 1, for Poisson noise ground truth was scaled to 
be between 0 and 5, and for speckle noise ground truth was scaled to be between 0 and 5 (with 3 looks 
in the noise generation). Eigure|^ shows reconstruction results for the baseline reconstruction and the 
best reparameterized reconstruction, and also lists the RMSE averaged over 20 noise realizations for each 
estimate; we see a significant improvement in terms RMSE and visual quality when reparameterizing. Eigure 
m shows the same reconstruction results on a zoomed-in view for more detail. 


7.2 Natural image experiments 

The Shepp-Logan experiments show value in reparameterizing in terms of RMSE and visual image quality. 
We now more extensively analyze reparameterizations through comparing reconstruction results on natural 
images to see which choice of reparameterization minimizes the empirical risk ( [2^ . The bank of images 
considered consists of n = 10 images from [7]. These are commonly seen in the image processing literature, 
and include: Barbara (512x512), Boat (512x512), Cameraman (256x256), Couple (512x512), Eingerprint 
(512x512), Hill (512x512), House (256x256), Lena (512x512), Man (512x512), and Peppers (256x256). 
There are m = 20 noisy realizations for each test image. Eor the Bernoulli case, the true images were scaled 
to have intensities on [0,1], and for the Poisson case, the images were scaled to have nonnegative intensities of 
mean 5 and 10 per pixel. Eor the speckle noise experiments, the true images were scaled to have nonnegative 
intensities of mean 3 and 5 per pixel, with S = 3 and 5 looks. 


7.2.1 Exponential family experiments 

Eigure shows plots of how the relative risk R{f,p) with respect to k for the Bernoulli and Poisson cases, 
using the proposed (15) and log-exp (18), respectively. These curved are averaged over the 10 images with 
20 noise realizations each and indicate optimal settings of the parameter k in each case. The data for these 
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plots can be found in Tables lAlllBTi and |Cl| in the appendices, where the averaged RMSEs for each image’s 
noise realizations are listed; the last row of the tables corresponds to the R plots in Figure We see that 
the baseline comparison performs worse than the reparameterizations considered. 

The appendices also show example reconstructions. Figure [AT] shows example reconstructions from one 
noise realization of the Hill and House data sets under Bernoulli noise. Figure |B1| shows example recon¬ 


structions from the Poisson realizations for the Barbara and Couple images with mean 5. Figure [CT] shows 
example reconstructions from the Poisson realizations for the Lena and Man images with mean 10. 

We note that for these noise examples, the baseline negative log-likelihood with respect to image intensities 
is convex, which fits well into the SpaRSA framework. 


7.2.2 Speckle noise experiments 

Figure]^ shows how R(/, p) changes as the reparameterization functions change with respect to k for speckle 
noise using the proposed ( [21] ). Averaged over the 10 images with 20 noise realizations each, we see that as 
k grows, R initially decreases but then stabilizes for large k. We remark that we search over small values 
of k (between 0 and 2) for practical reasons, as large k can cause numerical instabilities. Table D1 in the 


appendices lists the RMSE averages of the baseline comparison and k = 2 reparameterization for each image 
considered; we only show the values at /c = 2 because of the observed leveling off behavior. 

Figure [ dT] shows example reconstructions from one noise realization of speckle noise of the Peppers (S = 3) 
and Cameraman {S = 5) data sets with mean 3. Figure D2 shows example reconstructions from one noise 
realization of speckle noise of the Fingerprint {S = 3) and Boat {S = 5) data sets with mean 5. 

We note that the baseline negative log-likelihood data-fit for the speckle noise is not convex, but we can 
still apply SpaRSA to find a minimizer. 


7.2.3 Computational-statistical tradeoff 

We briefly revisit the idea of a computational-statistical tradeoff discussed in SectionFigure [^illustrates 
one example from Poisson noise (specifically, we are looking at results from the Couple image with mean 
10). The plots show how RMSE and the number of SpaRSA iterations behave with the parameter k. We see 
that RMSE vs. k has a dip and that Iterations vs. k has a generally increasing behavior. Plotting RMSE 
vs. iterations shows us how well we can expect to perform given a limited computational budget. 


7.3 Discussion 


In these experiments, the values of r were chosen clarivoyantly to minimize the RMSE for each reparame¬ 
terization considered. In practical settings, we do not have access to ground truth to compute RMSE, but 
we can use cross-validation to select r and the specific k parameter within the family considered to obtain 
an estimate. If on a time budget, we can leverage the experimental results to eliminate the search over Us: 


• In the Bernoulli case. Table |A1| in the appendices shows that the value of k that minimizes the risk 
R(/, p) for the average of the bank of 10 images is k = 0.3. 


In the Poisson case. Tables |B1| and |C1| in the appendices show that k = 2 minimizes the average risk 
R(/, p) over the 10 images with mean 5, and k = I minimizes the average risk for images with mean 10; 
to express /c as a fraction of the mean we can say k = 0.4x*. In practical settings, for a large number 
of pixels we can roughly estimate k ~ 0.4^ to use in the data-fit (18). 


In the speckle noise case, we have already established that as k grows, the average risk R(/, p) stabilizes. 
In practical settings, we choose to use k = 2 when applying it to the data-fit (21). 


8 Conclusion 

We have introduced a reparameterization framework for three types of noise commonly found in imaging 
problems, where the reparameterizations proposed are determined by the structure of the negative log- 
likelihood. This essentially models the image as a parameterized function, allowing us to implicitly modify the 
regularizer. Image reconstruction is still practical with this change, as we can still use well-studied proximal 
denoising subroutines in the optimization. Through examples from the exponential family and speckle noise. 
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we see that reparameterization yields gains in RMSE, and has computational implications. Because searching 
over all possible reparameterization functions is not practical, the functions presented in this paper are by 
no means “optimal”; we simply introduce a generalized framework (with example reparameterizations) that 
may be improved upon with the selection of other functions in future work. 
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(a) Truth (b) Noisy observations (c) Baseline recovery (d) Reparameterized 

(Bernoulli) (RMSE: 0.1554) recovery, k = 0.05 

(RMSE: 0.1249) 



(e) Truth (f) Noisy observations (g) Baseline recovery (h) Reparameterized re- 

(Poisson) (RMSE: 0.2074) covery, A; = 4.0 (RMSE: 

0.1442) 



(i) Truth (j) Noisy observations (k) Baseline recovery (1) Reparameterized re- 

(Speckle) (RMSE: 0.6483) covery, k = 2.0 (RMSE: 

0.3229) 


Figure 4: Shepp-Logan phantom experiments: Denoising results using Bernoulli noise (first row), Poisson noise 
(second row), and speckle noise (third row). We show the ground truth, noisy observations, the baseline recovery 
(24), and the best reparameterized recovery. The RMSEs of our reconstructions (averaged over 20 noise realizations) 
are also listed. 
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(a) Truth (b) Noisy observations (c) Baseline recovery 

(Bernoulli) (RMSE: 0.1554) 



(d) Reparameterized 
recovery, k = 0.05 

(RMSE: 0.1249) 



(e) Truth 


(f) Noisy observations (g) Baseline recovery (h) Reparameterized re- 
(Poisson) (RMSE: 0.2074) covery, A; = 4.0 (RMSE: 

0.1442) 



(i) Truth (j) Noisy observations (k) Baseline recovery (1) Reparameterized re- 

(Speckle) (RMSE: 0.6483) covery, k = 2.0 (RMSE: 

0.3229) 


Figure 5: Shepp-Logan phantom experiments (zoomed-in view): Denoising results using Bernoulli noise (first 
row), Poisson noise (second row), and speckle noise (third row). We show the ground truth, noisy observations, the 
baseline recovery and the best reparameterized recovery. The RMSEs of our reconstructions (averaged over 20 
noise realizations) are also listed; note that the RMSEs are for the entire image and not for the zoomed-in parts. 
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0.16 


0.155 


^ 0.15 

0.145 


O.I4L 

0 



■ ■■■Baseline 
— Reparameterization 


0.1 0.2 0.3 0.4 0.5 

k 

(a) Bernoulli 



k 


(b) Poisson (mean 5) 



k 

(c) Poisson (mean 10) 


Figure 6: Bernoulli and Poisson noise i^(/, p) plots: These plots show how the relative risk R{f^ p) for 
the proposed reparameterizations (solid lines) varies as the parameter k changes for Bernoulli and Poisson 
noise, using the proposed (15) and log-exp ( [l8| ), respectively. Also shown are ^(/,p)’s for the baseline 
reconstructions 


2^ (dashed lines) as comparisons. These plots are averages over 10 images with 20 noise 
realizations each; the data used to generate these plots can be found in Tables |A1[ |B1[ and |C1| in the 
appendices. 




)) 



■ ■■■S=3, Baseline 
— S=3, Reparam. 

■ ■■■S=5, Baseline 
— S=5, Reparam. 


(a) Mean 3 


(b) Mean 5 


Figure 7: Speckle noise R{f^p) plots: These plots show how the relative risk R{f^p) for the proposed 
reparameterization (solid lines) varies as the parameter k changes for various combinations of mean image 
intensity and number of looks S', using the proposed ( [^. We see that R{f, p) levels off as k increases. 
Also shown are ^(/, p)’s for the baseline reconstructionsT2^ (dashed lines) as comparisons. These plots are 
averages over 10 images with 20 noise realizations each. Table [DT] in the appendices shows the RMSE values 
for the baseline comparison and the k = 2 reconstruction for each image considered. 
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35 40 

Iterations 
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50 


(a) RMSE vs. k 


(b) Iterations vs. k 


(c) RMSE vs. Iterations 


Figure 8: Computational-statistical tradeoff: These plots show algorithmic performances for the Couple 
image with mean 10 corrupted by Poisson noise using the proposed log-exp ( p!8| ). A dip in the RMSE vs. k 
plot (a) and a generally increasing behavior in algorithm Iterations vs. k (b) shows us that there is a dip in 
RMSE vs. Iterations (c). This plot in (c) suggests how well (in terms of RMSE) we can expect to perform 
with the proposed reparameterizations when given a computational budget. 
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A Additional numerical results for Bernoulli experiments 


Image 

Baseline 

0.05 

0.10 

0.15 

0.20 

k 

0.25 

0.30 

0.35 

0.40 

0.45 

0.50 

Barbara 

0.1766 

0.1685 

0.1661 

0.1643 

0.1635 

0.1636 

0.1632 

0.1632 

0.1634 

0.1635 

0.1635 

Boat 

0.1507 

0.1415 

0.1391 

0.1373 

0.1368 

0.1375 

0.1372 

0.1379 

0.1376 

0.1379 

0.1376 

Cameraman 

0.1658 

0.1609 

0.1585 

0.1573 

0.1581 

0.1590 

0.1600 

0.1605 

0.1610 

0.1613 

0.1614 

Couple 

0.1613 

0.1531 

0.1509 

0.1498 

0.1494 

0.1503 

0.1503 

0.1508 

0.1507 

0.1509 

0.1509 

Fingerprint 

0.2207 

0.2111 

0.2087 

0.2056 

0.2029 

0.2010 

0.1998 

0.1995 

0.1990 

0.1991 

0.1991 

Hill 

0.1443 

0.1373 

0.1351 

0.1331 

0.1322 

0.1316 

0.1314 

0.1313 

0.1313 

0.1314 

0.1315 

House 

0.1148 

0.1099 

0.1092 

0.1104 

0.1104 

0.1095 

0.1102 

0.1103 

0.1107 

0.1111 

0.1112 

Lena 

0.1254 

0.1165 

0.1153 

0.1145 

0.1141 

0.1159 

0.1147 

0.1146 

0.1151 

0.1152 

0.1158 

Man 

0.1542 

0.1444 

0.1437 

0.1424 

0.1412 

0.1412 

0.1411 

0.1404 

0.1406 

0.1408 

0.1409 

Peppers 

0.1754 

0.1678 

0.1651 

0.1638 

0.1628 

0.1627 

0.1630 

0.1628 

0.1626 

0.1629 

0.1631 

Average 

0.1589 

0.1511 

0.1492 

0.1479 

0.1471 

0.1472 

0.1471 

0.1471 

0.1472 

0.1474 

0.1475 


Table Al: Bernoulli noise data: RMSE results (all averaged over 20 noise realizations) for each image for the 
baseline reconstruction (24) and for 10 different reparameterizations from k = 0.05 to 0.50 with the proposed ( p^ . 
Each row has the best performing RMSE in bold. The last row averages the columns, which is the relative risk R{f, p) 
defined by 0 . and these values are plotted in Eigure j^ We see that the baseline estimate performs worse than the 
best reparameterization (in fact, for all the k values chosen, the baseline is overall worse). 




(a) Truth 


(b) Observations 



(c) Baseline recovery (d) Reparameterized re¬ 
covery {k = 0.40) 



(e) Truth 


(f) Observations 



(g) Baseline recovery (h) Reparameterized re¬ 
covery {k = 0.10) 



Figure Al: Bernoulli noise images: Reconstruction results from the “Hill” image (top row); reconstruction 
results from the “House” image (bottom row). These images are all displayed on the same intensity scale (from 0 to 
1.1 X max(x*)). 
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B Additional numerical results for Poisson experiments (mean 5) 


Image 

Baseline 

0.5 

1.0 

1.5 

2.0 

k 

2.5 

3.0 

3.5 

4.0 

4.5 

5.0 

Barbara 

0.1472 

0.1404 

0.1385 

0.1375 

0.1369 

0.1365 

0.1365 

0.1367 

0.1369 

0.1369 

0.1370 

Boat 

0.1091 

0.1045 

0.1030 

0.1029 

0.1030 

0.1035 

0.1038 

0.1042 

0.1043 

0.1046 

0.1048 

Cameraman 

0.1228 

0.1115 

0.1107 

0.1105 

0.1109 

0.1114 

0.1120 

0.1129 

0.1134 

0.1140 

0.1144 

Couple 

0.1171 

0.1130 

0.1119 

0.1121 

0.1116 

0.1120 

0.1123 

0.1126 

0.1129 

0.1134 

0.1136 

Fingerprint 

0.1698 

0.1599 

0.1557 

0.1528 

0.1515 

0.1513 

0.1516 

0.1521 

0.1520 

0.1525 

0.1527 

Hill 

0.1110 

0.1039 

0.1019 

0.1012 

0.1009 

0.1010 

0.1011 

0.1015 

0.1016 

0.1018 

0.1022 

House 

0.0931 

0.0888 

0.0868 

0.0854 

0.0854 

0.0852 

0.0856 

0.0856 

0.0859 

0.0860 

0.0862 

Lena 

0.0889 

0.0877 

0.0850 

0.0838 

0.0838 

0.0840 

0.0846 

0.0849 

0.0852 

0.0856 

0.0860 

Man 

0.1154 

0.1091 

0.1066 

0.1057 

0.1050 

0.1050 

0.1053 

0.1054 

0.1054 

0.1056 

0.1060 

Peppers 

0.1131 

0.1061 

0.1046 

0.1044 

0.1044 

0.1050 

0.1055 

0.1059 

0.1063 

0.1068 

0.1074 

Average 

0.1187 

0.1125 

0.1105 

0.1096 

0.1093 

0.1095 

0.1098 

0.1102 

0.1104 

0.1107 

0.1110 


Table Bl: Poisson noise data (mean 5): RMSE results (all averaged over 20 noise realizations) for each image for 
the baseline reconstruction |24| and for 10 different reparameterizations from /c = 0.5 to 5.0 with the proposed log-exp 


T^. Each row has the best performing RMSE in bold. The l ast row averages the columns, which is the relative risk 


R{f,p) defined by 0- and these values are plotted in Eigure|^ We see that the baseline estimate performs worse 


than the best reparameterization (in fact, for all the k values chosen, the baseline is overall worse) 



(a) Truth (b) Observations (c) Baseline recovery (d) Reparameterized re¬ 

covery {k = 3.0) 



(e) Truth (f) Observations (g) Baseline recovery (h) Reparameterized re¬ 

covery (k = 2.0) 


Figure Bl: Poisson noise images (mean 5): Reconstruction results from the “Barbara” image (top row); recon¬ 
struction results from the (zoomed-in) “Couple” image (bottom row). These images are all displayed on the same 
intensity scale (from 0 to 1.1 x max(x*)). 
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C Additional numerical results for Poisson experiments (mean 10) 


Image 

Baseline 

1 

2 

3 

4 

k 

5 

6 

7 

8 

9 

10 

Barbara 

0.1369 

0.1271 

0.1246 

0.1228 

0.1221 

0.1220 

0.1220 

0.1221 

0.1222 

0.1223 

0.1224 

Boat 

0.0954 

0.0906 

0.0890 

0.0889 

0.0888 

0.0888 

0.0891 

0.0893 

0.0895 

0.0897 

0.0900 

Cameraman 

0.1022 

0.0955 

0.0924 

0.0924 

0.0928 

0.0933 

0.0938 

0.0942 

0.0947 

0.0951 

0.0955 

Couple 

0.1031 

0.0985 

0.0966 

0.0965 

0.0964 

0.0966 

0.0968 

0.0971 

0.0973 

0.0975 

0.0978 

Fingerprint 

0.1347 

0.1301 

0.1280 

0.1263 

0.1255 

0.1255 

0.1257 

0.1259 

0.1260 

0.1262 

0.1264 

Hill 

0.0941 

0.0905 

0.0895 

0.0887 

0.0884 

0.0882 

0.0884 

0.0885 

0.0887 

0.0889 

0.0891 

House 

0.0745 

0.0717 

0.0711 

0.0709 

0.0708 

0.0704 

0.0705 

0.0706 

0.0707 

0.0708 

0.0709 

Lena 

0.0755 

0.0727 

0.0717 

0.0716 

0.0717 

0.0718 

0.0721 

0.0723 

0.0726 

0.0729 

0.0731 

Man 

0.0975 

0.0943 

0.0925 

0.0916 

0.0910 

0.0908 

0.0909 

0.0910 

0.0911 

0.0913 

0.0914 

Peppers 

0.0937 

0.0878 

0.0873 

0.0865 

0.0865 

0.0866 

0.0869 

0.0872 

0.0875 

0.0878 

0.0882 

Average 

0.1008 

0.0959 

0.0943 

0.0936 

0.0934 

0.0934 

0.0936 

0.0938 

0.0940 

0.0942 

0.0945 


Table Cl: Poisson noise data (mean 10): RMSE results (all averaged over 20 noise realizations) for each image 
for the baseline reconstruction \2A\ and for 10 different reparameterizations from /c = 1 to 10 with the proposed 
log-exp ( |18| ). Each row has the best performing RMSE in bold. The last row averages the columns, which is the 
relative risk R{f,p) dehned by 0- and these values are plotted in Eigure We see that the baseline estimate 
performs worse than the best reparameterization (in fact, for all the k values chosen, the baseline is overall worse). 



(a) Truth (b) Observations (c) Baseline recovery (d) Reparameterized re¬ 

covery (k = 3) 



(e) Truth (f) Observations (g) Baseline recovery (h) Reparameterized re¬ 

covery (k = 5) 


Figure Cl: Poisson noise images (mean 10): Reconstruction results from the “Lena” image (top row); recon¬ 
struction results from the “Man” image (bottom row). These images are all displayed on the same intensity scale 
(from 0 to 1.1 X max(x*)). 
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D Numerical results for speckle noise experiments 


Image 

Mean 3, 
Baseline 

S = 3 
k = 2 

Mean 3, 
Baseline 

S = 5 
k = 2 

Mean 5, 
Baseline 

S = 3 
k = 2 

Mean 5, 
Baseline 

S = 5 
k = 2 

Barbara 

0.1717 

0.1491 

0.1664 

0.1397 

0.1969 

0.1493 

0.1810 

0.1399 

Boat 

0.1373 

0.1183 

0.1221 

0.1056 

0.1568 

0.1184 

0.1385 

0.1059 

Cameraman 

0.1530 

0.1312 

0.1398 

0.1155 

0.1716 

0.1313 

0.1594 

0.1161 

Couple 

0.2176 

0.1273 

0.1694 

0.1140 

0.1794 

0.1271 

0.1656 

0.1139 

Fingerprint 

0.2261 

0.1862 

0.1980 

0.1607 

0.2595 

0.1861 

0.2120 

0.1606 

Hill 

0.1580 

0.1157 

0.1446 

0.1047 

0.1649 

0.1160 

0.1517 

0.1049 

House 

0.1504 

0.1023 

0.1362 

0.0896 

0.1441 

0.1012 

0.1369 

0.0891 

Lena 

0.1331 

0.0981 

0.1219 

0.0870 

0.1413 

0.0982 

0.1291 

0.0868 

Man 

0.1530 

0.1218 

0.1412 

0.1092 

0.1656 

0.1218 

0.1543 

0.1093 

Peppers 

0.1489 

0.1256 

0.1329 

0.1084 

0.1764 

0.1247 

0.1629 

0.1077 

Average 

0.1649 

0.1276 

0.1472 

0.1134 

0.1757 

0.1274 

0.1591 

0.1134 


Table Dl: Speckle noise data: RMSE results (all averaged over 20 noise realizations) for each image for the 
baseline reconstruction and for the reparameterization with k — 2 with the proposed ( [^ ; we only state the 
results ioT k — 2 because Figure [^suggests that the RMSE levels off as k increases (and for the sake of space). The 
last row averages the columns, which is the relative risk R(/, p) defined by 0. and these values are plotted in Figure 
We see that the baseline estimate performs worse than the best reparameterization at the leveled-off valued of 
k = 2. 



(a) Truth (b) Observations (c) Baseline recovery (d) Reparameterized re¬ 

covery {k = 2) 



(e) Truth (f) Observations (g) Baseline recovery (h) Reparameterized re¬ 

covery (k = 2) 


Figure Dl: Speckle noise images (mean 3): Reconstruction results from the “Peppers” image with S = 3 (top 
row); reconstruction results from the “Cameraman” image with S' = 5 (bottom row). These images are all displayed 
on the same intensity scale (from 0 to 1.1 x max(x*)). 
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(a) Truth (b) Observations (c) Baseline recovery (d) Reparameterized re¬ 

covery (k = 2) 



(e) Truth (f) Observations (g) Baseline recovery (h) Reparameterized re¬ 

covery (k = 2) 


Figure D2: Speckle noise images (mean 5): Reconstruction results from the “Fingerprint” image with S = 3 
(top row); reconstruction results from the “Boat” image with S' = 5 (bottom row). These images are all displayed on 
the same intensity scale (from 0 to 1.1 x max(x*)). 
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